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SUMMARY 


An accurate and efficient numerical solution algorithm is established 
for solution of the high Reynolds number limit of the Navier-Stokes equations 
governing the multi -dimensional flow of a compressible essentially inviscid 
fluid. The theoretical basis employs finite element interpolation theory 
within a dissipative, formulation established using Galerkin criteria within 
the Method of Weighted Residuals (MWR). An implicit iterative solution al- 
gorithm is developed, employing tensor product bases within a fractional 
steps integration procedure, that significantly enhances solution economy 
concurrent with sharply reduced computer hardware demands. The algorithm is 
evaluated for resolution of steep field gradients and coarse grid accuracy 
using both linear and quadratic tensor product interpolation bases. Numeri- 
cal solutions for linear and non-linear, one-, two- and three-dimensional 
examples confirm and extend the linearized theoretical analyses, and results 
are compared to competitive finite differenced-derived algorithms. 



INTRODUCTION 


Finite element concepts burst upon the computational fluid dynamics 
scene about a decade ago in the guise of a triangle. The primary motiva- 
tion was the profuse'' geometric flexibility, in contrast to the then- 
current finite difference limitation to ‘regular grids. In the ensuing 
interval, the inherent versatility of the basic finite element concept has 
proven difficult to master v/ithin a computationally economical framework. 

In the same period, development of regularizing coordinate transformation 
(ref. 1, 2) has markedly extended the applicability of efficient finite 

/ 

difference recursion formulae to non-regular shaped solution domain closures. 

What was patently obscure in the early work, but is becoming convin- 
cingly transparent, is that the finite element/ weighted residuals theoretical 
basis provides a foundation for derivation of optimal ly-accurate (in the 
appropriate norm) numerical algorithms for solution of general categories in 
fluid mechanics. The theoretical support for finite element solution of 
linear elliptic equations is complete, and in particular one is assured that 
a finite element potential flow solution is optimally accurate in the 
(energy) norm in comparison to all other methods (ref, 3). Solution of ini- 
tial-valued problem- descriptions is quite typical in fluid mechanics, and in 
the L 2 norm the finite element algorithm is confirmed optimally accurate for 
linear parabolic equations (ref^ 4). The numerical extension to non-linear 
parabolic equations, as appropriate for boundary layer flows, has confirmed 
extremization of the energy norm for both laminar (ref. 5) and turbulent 
flows (ref. 6). The latter is of particular interest, since the energy norm 
is a strongly non-linear function of the mean flow field gradients through the 
effective turbulent viscosity. While by no means constituting a theoretical 
proof, these results to indicate that the linear finite element theory may be 
extensible to the more interesting problem classes. 
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A particularly difficult problem class in fluid dynamics corresponds 
to high Reynolds number flow prediction. The governing Navier-Stokes 
equations are generally a non-linear elliptic boundary value description, 
but the importance of the viscosity term is com.pletely dwarfed by the non- 
linear convective acceleration everywhere away from a wall. Furthermore, 
the continuity equation possesses no viscosity-like term, hence exhibits 
uniformly the hyperbolic description pervading the entire flowfield region 
away from walls. The primary objective of the research reported herein is 
to derive and numerically evaluate a finite element/weighted residuals solu- 
tion algorithm applicable to large Reynolds number flow prediction. It was 
decided that the algorithm must be implicit so as to handle physical viscosity 
effects as appropriate. However, primary emphasis rests on determination of 
accuracy and convergence phenomena for dominantly inviscid forms of the 
Navier-Stokes equations, ie. the Euler equations. Since three-dimensional 
flow field prediction is the eventual goal, the developed algorithm must be 
computer core and CPU efficient. Furthermore, since multi-dimensional fluid 
mechanics predictions are of necessity almost universally performed on coarse 
computational grids, particular emphasis is placed on coarse-grid, accuracy 
assessment. The accurate and efficient finite element tensor-product algorithm 
that has been derived to meet these requirements is reported herein. 

PROBLEM STATEMENT 

The basic requirement is to establish a numerical solution algorithm to 
accurately and efficiently model the substantial time derivative associated 
with all flow field descriptions (save potential flow). This hyperbolic operator 
dominates high Reynolds number flows, as governed by the Navier-Stokes equations. 
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The conservation form of this familiar partial differential equation system is 


L(p) = If + si- ^PU-;) - 0 
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In equations (l)-(3), p is the density, pu^ the momentum vector, e is the 
specific energy, and p is the static pressure defined by the equation of 
state. 


P = 



■rr pu.u. 
2 ^11 


(4) 


where y is the ratio of specific heats. For present purposes, the heat flux 
vector 8- and the Stokes stress tensor cr. . are both assumed negligible; hence, 

o I J 

equations (l)-(3) are hyperbolic. Equations (2), l<i£3, are explicitly non- 
linear, while equations {!) and (3) are quasi-linear in the expressed dependent 
variable, and each require characterization. 

Equations (l)-(3) describe the transient evolution of the element of 
PP n-dimensional space spanned by the x^ coordinate system, 
l<i<n. The domain of the solution is e x t e x- xOt^jt) with closure 
9fi E R*^ ^ X t. On n, each member of {q} is the solution to 


= It " h'l * { 


= 0 


(5) 


where f(q|,) is specified in equations (2)-{3). All applicable boundary con- 
ditions on 9fl are contained within the expression 

it(q) = aiq + 92 Ij- n^ + 33 = 0 


( 6 ) 



where n- is the local outward-pointing unit normal vector and the a. are 
J T 

specified coefficients. An initial condition on = r” x t is required, 

0 0 

hence 

q(x^., t^) = qQ(x.) (7) , 

'Since regularizing coordinate transformations are available, no generality 
is lost in assuming x^. a Cartesian coordinate system spanning R^. 

FINITE ELEMENT SOLUTION ALGORITHM 

A dissipative finite element solution algorithm is established for 
equations (5)-(7), hence the Navier-Stokes and Euler equations. Assuming 
Ufig = n, the domain of L(q), that the are non-overlapping, and that 
^ ^ where UR^ is the finite element discretization of R*^, let each 

member of {q} be interpolated on as ‘ 

qg{x., t) = {Nj,(x^.)}^{Q(t)}g (8) 

The elements of {Nj^(x^-)} are polynomials in x^, complete to degree k, and 
form a cardinal basis (ref. 7). The expansion coefficients {Q(t)}g are un- 
known, a solution algorithm for which is required established to determine the 
temporal evolution of the dependent variable system {q}. To accomplish this, 
substitute equation (8) into equations (5)-(6), and set to zero the integral 
of each over R^ and SR^ after weighting by {N|^(x^. )}. In addition, as sug- 
gested in reference 8,. set the weighted integral of the vector gradient of 
equation (5) to zero. Identifying the vector and scalar multipliers 8^ and' X, 
which must be determined, combine these expressions into the matrix equation 
system. 
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Sg f {N}L(qg)dT - ej {N} ^L{q JdT + xf {NH(q )da ={0} (9) 

Re R? SR^naR" 

Equations (9) are systems. of ordinary differential equations written on the 

temporal evolution of the discrete approximation Z{Q}g to each qg(x^- , t) of 

e 

{q}. These differential equations are uncoupled in the temporal .derivatives , 
and Sg is the familiar finite element assembly operator mapping local opera- 
tions to the global reference frame (cf., ref, 4). The scalar multiplier X 
is conventionally employed to enforce the discretized boundary condition state- 
ment, equation (6). 

The B.J are scalar components of an n-dimensional vector, the determina- 
tion of which is required. This term represents the additional requirement 
that the gradient of the solution error in L(qg) be orthogonal to the inter- 
polation basis {Nj,}. The desired form for the augmented MWR algorithm is 
achieved using a Green-Gauss theorem. Letting p. = v-A , where A is the 
measure of the finite element domain R^, yields 


' 

•'nn 


1 + A, 


V. {N} L(q )dT + A {N}£(q )da = {0} 

^ J 9R^H3R" 


The n scalar components of v. can be element-dependent in the general case, 

«J 

and the closure surface integral stemming from use of the Green-Gauss theorem 

vanishes identically. Independent of the dimension n of R*^, equation (10) 

yields an ordinary differential equation system for solution of 2{Q(t)} = {Q}, 

e ® 

of the form 





The superscript prime denotes (ordinary) differentiation with respect to 
time. The first two' terms, common for all qj^, see equations (5)-(6), are 
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(13) 


Here, {U^} is the nodal distribution of the discrete representation of the 
velocity field u^(x,t). The source term {f(q)}g is distinct for each qj^, 
and represents a non-homogeneity at the minimum. 

As mentioned, an implicit integration algorithm for equation (11) is re- 
quired. A familiar single-step procedure is 



(14) 


where j is the time-step ’index, h is the integration stepsize, and 0 is a 
parameter 0 < e <_ 1 controlling implicitness. Following the usual manipula- 
tions (ref. 4), insertion of equation (11) into (14) yields a large order, non- 
linear algebraic equation system. The Newton matrix iteration algorithm for 
solution of this system is 





The dependent variable is the iteration vector, 


(15) 




( 16 ) 



where p is the iteration index. The right side of equation (15) is the 


4* W 

homogeneous form of equation (14) evaluated with the p~^ iterate 






( 17 ) 


where 


{Se>? = Me«}? + (f> 


e 


( 18 ) 


Note that equations (17)-(18) are defined solely in terms of inner products on 
elements, with the assembly operator yielding the equivalent global contri- 
bution. The vanishing of {F} to within definition of a computed zero yields 
equation (15) homogeneous, hence convergence of the iteration process. By 
definition, the Jacobian is the derivative of equation (17) with respect to 
Hence, 


[J] = s, 


[cjp + h0[u]^ + he 


Jm 


+ {f} 
e e 


mF: 


(Q) 




(19) 


where the final term accounts for contributions stemming from explicit non- 
linearity. ATI operations involve matrix inner products of an elemental basis, 
hence implicitly independent of the dimension n of R*^. The rank of [J] at 
least equals the order of {6Q}; specific (Dirichlet) boundary constraints are 
applied within the evaluation of {F}. 

As opposed to the conventional use of multidimensional, finite element 
interpolation functions {N(x^)}, the three-dimensional requirements demand 
a spatial factorization that permits replacement of the large, sparse-matrix 
Jacobian operations with elementary banded-matrix procedures. The theoreti- 
cal operations are to replace the multi -dimensional interpolation bases with 
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tensor product bases (ref. 9), and to Implement the method of fractional 
steps (ref. 10) to resolve solution operations onto scalar components 
parallel to coordinate axis. Hence, the interpolation basis in equation 
(8) becomes, for three-dimensional space 

{N^(x^ )} ^ ® ® (20) 

where @ signifies the tensor product. The equivalent tensor matrix products 
are similarly expressed; for example, the matrix equivalent of the initial - 
value operator, equation (12) becomes 

Me <^1) 

where 

ty e " L J[l " ''e '’a) 

e *- -* ■ 

The similar operations for the convection operator, equation (13) yields 

Me = Lit " '^e ^ ''.]«k>| sLi 'yI«kJ'"k)L^a 

L a J aL 

and a is not a tensor summation index in equations (22)- (23). 

With the finite element matrix equivalents of the terms in equation (11) 
recast as tensor products, the method of fractional steps (ref, 10) is employed 
to establish the desired operations for the Newton matrix iteration algorithm, 
equation (15). It is elementary operation to evaluate the Jacobian as the 
tensor matrix product, 

[J]=^ [JJ0[O,]0[Jb] - (24) 
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and each of the [J^] exhibit the desired banded matrix structure. Spe- 
cifically, for k = 1 in equation (20), is tridiagonal, while for k = 2 
it is dominantly pentadi agonal . Hence, the tensor matrix solution algorithm 
equation (15) becomes 

[j=«3)5+i] ® = 

- [f3(Q)?+J C26) 

Here, Q and Q represent intermediate iterates and 6Q is interpreted as the 
iteration vector for each respective iterate. 

THEORETICAL ANALYSIS 

A von Neumann stability analysis (ref. 11) can quantize the formal 
order of accuracy of the develoepd algorithm for a linearized one-dimensional 
equation, hence predict an appropriate value for v. Therefore, consider the 
x- momentum equation with constant advection velocity U^, ie 

The analytical solution to equation (26) is the Fourier expansion 

u(x,t) = V exp w(x - U^t)J (27) 

where i e /T, w = 2frA is the wave number where X is wavelength and V is 
the initial velocity distribution u(x,0.). This solution corresponds to the 
diffusion- and dispersion-free advection of the initial wave form parallel to 
the x*^ axis with velocity 
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The Fourier analysis of the discrete algorithm is best accomplished 
using equation (11). It is- readily facilitated only for specification of 
k = 1 in equation (20), wherein the assembly operation yields an elementary 
recursion relation form of the algorithm. The resultant expression for 
equation (22) is 



'2 r 

< 

r-i -ii 


_1 2_ 

4* ^ 

I 

h~» 

ilrL- 


(28) 


for V a constant. For a uniform discretization, the assembly of the first 
term of equation (11) over the two elements sharing node j yields. 




a 


(l+3v) + 4Q^ + (l-3v)Qj^j 


(29) 


where is the uniform measure of the discretization in the direction of x^. 
Similarly, the convection term in equation (23) becomes 

For the constant advecti on. velocity,' {11'^}^ = 

The second form of equation (31) emphasizes the action of the dissipative al- 
gorithm in introduction of a viscosity-like term, ie. the difference operation 




“o 


-(l+2v)Qj_j + 4vQj + 


^T\ 




+ U^v 


■‘* 3-1 ^**3 ■ ‘* 3+1 


(31a) 


(31b) 
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equivalent of a second derivative term in equation (26) except for the omis- 
sion of a power of the measure Alternatively, the first form suggests 
the action of v 0 as a sort of upwind difference operation. It is import- 
ant to note, however, that these elementary interpretation are valid only 

for the selected case, ie. a constant and uniform. 

0 e 

The von Neumann stability analysis assumes the solution for the semi -dis- 
crete equivalent of the continuous Fourier form as 


u*(Ax,t) = Z Ug 



(32) 


where Ax =-a“, the uniform discretization of r“, j is the node indicator, and 
X = B + 15, where B and '5 are real numbers. Comparing equations (27) and 
(32), a difference between 3 and constitutes a disparity in the phase 
speed of propagation of V, hence phase error in .the discrete solution. Cor- 
respondingly, 6 0 introduces a real exponential argument yielding a damping 

(or growth) of the amplitude of the initial distribution V. Direct substi- 
tution of equations (29) and (31) into (32) and expanding the resultant ex- 
pressions for 3 and 5 in a Taylor series yields (ref. 12) 


3 = U 


1 + 




1 + V 

'180 


0(d«) 


(33) 


6 = 


vd 

■ ir* 



(34) 


Here, d = toAx and 0( ) indicates the order of the truncated term. Since 3 is 
the real component of X, it can be made identical to to order (Ax)® by re- 
quiring V ^ = /IF. Then, the phase accuracy of the discrete solution u* 
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agrees with to sixth order accuracy In Ax, i.e,, 0(Ax)®. For v 0, 

6 < 0 and an artificial damping is introduced. The resultant specific form 
for equation (32) is . 


r r 


u*(Ax,t) = V exp 


1 0 ) 


jAx - + 0(Ax®)j 


exp 


- 0 ) 


'kt + 0 (Ax' 


>] 


(35) 


where k e U^(Ax)V12v^ is the damping coefficient. Note that the damping is 
quite selective, occuring only for sufficiently large wave numbers w (small 
wave lengths) due to the u** factor. 

Two additional comments are warranted. Setting v = 0 eliminates 6 and 3 
is a fourth order accurate representation of the differential equation. This 
high order accuracy accrues with use of the simplest linear interpolation. 

It is obvious that the convection term in equation (31b) is the central dif- 
ference equivalent; the improvement to fourth order results directly from the 
finite element derived form for [C ie. equation (29). The normalized 
(1,4,1) weighting on the derivatives corresponds identically with a spline 
interpolation. The conventional finite difference practice, eg, reference 13, 
Is to replace equation (19) with a'^Q. which yields directly a degradation to 
overall second order accuracy. Secondly, the conventional finite difference 
practice to introduce dissipation is to add the "artificial viscosity" term 
p 3^u/8x^ to the parent differential equation (29). Repeating the semi-dis- 
crete Fourier analysis with this added term, and setting v e 0 but retaining 
the finite element derived initial-value term, yields (ref. 12). 

u*(Ax,t) = V exp 


1 w 


jAx - 


Uq + 0(Ax'*) 


expjj co^yt + 0(Ax^)J (36) 
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The basic fourth order phase accuracy remains intact, but the artificial 
damping is less selective due to the diminution of the wave number ex- 
ponent to 2 and the appearance of the term of order Ax^. 

The stability of the algorithm for the linearized equation can be as- 
sessed using a fully discrete Fourier analysis of equation (14) combined 
with (11). In this instance, letting n denote the time index, the approximate 
solution form is 


u**^(Ax,At) = g*^exp QwjAx^ (37) 

The form of the amplification matrix g*^ is sought, since the discrete solution 
will propagate and damp/grow dependent upon its real and imaginary arguments. 
Previous analyses (ref. 14), for the non-dissipative algorithm, indicate that 
the trapezoidal rule (Q e %) is the sole suitable selection. Retaining this 
definition, the amplification factor for the dissipative finite element al- 
gorithm is 

l+% cosCwAx) - 3Cv sin^(^Ax) - i|-(C+2v)sin(coAx) 

g = 4 (3g) 

l+% cos(tjoAx) + 3Cv sin^(%oAx) + i^C-2v)sin(uAx) 

C E U^Ax/At is the Courant Number, and the numerator and denominator of g are 
complex conjugates for v e 0. Therefore, the basic non-dissipative algorithm 
is neutrally stable, ie. jg] = 1 for all Ax and At, hence error induced by 
the solution will propagate undamped and unmagnified throughout the solution 
domain. Selecting v 0 destroys this neutral stability; therefore, define g 
in terms of real and imaginary parts as 

g = Y + ir (39) 
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Hence, y > 0 quantizes the dissipative mechanism, while r determines the 
phase accuracy. In particular. 



yields the normalized phase velocity of the approximate solution u**^. Figure 
1 graphs equation (40) for the non-dissipative form; the analytical solution 
corresponds to the horizontal line at u*/U^ = 1. The solid curves represent 
the finite element algorithm, while the dashed curves correspond to equation 
(38) as modified by the finite difference diagonalized initial-value matrix 
form. The superior performance of the present theory is clearly evident. At 
a modest Courant number, eg. 0.01 < C < 0.5, the algorithm accurately resolves 
all wavelengths X > 5Ax, while for the latter this occurs' only for X > 15Ax. 

The action of the added dissipative 'mechanism (v ^0) is to improve the 
phase accuracy in the short wave length region. Table 1 lists select evalua- 
tions of equation (38) for C = 0.25 forvarious v and X. The dissipation 
level pis defined as 


p = ^ Inlgj (41) 

It is evident that v exerts a profound correction to 9 on the interval 
2 <^X £5, hence produces a closer approximation to the correct solution re- 
garding phase accuracy. The penalty for improved phase accuracy is the 
corresponding introduction of dissipation. The dissipation level p is modestly 
sensitive to X and nearly linearly dependent on v. 
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PHASE CELERITY - U*/U 




Influence of Dissipation Level v on Phase Celerity and Dissipation Level, C = 0.25 


Wavelength 

X(nAx) 


Phase Celerity = 0 



Dissipation Level - 

U 

v=0. 

v=0.1 

v=0.25 

.v=0.5 

v=0 

v=0.1 

v=0.25 

v=0.5 

2.00 

0.00000 

0.00000 

0.00000 

0.00001 

0.0 

0.01241 

0.03438 

0.07991 

2.25 

0.34486 

0.38455 

0.59418 

1.26043 

0.0 

0.01615 

0.03854 

0.06108 

2.50 

0.58252 

0.62651 

0.82294 

1.26599 

0.0 

0.01637 

0.03611 

0.04907 

2.75 

0.72707 

- 0.76034 

0.90723 

1.20781 

0.0 

0.01547 

0.03322 

0.04358 

3.00 

0.81442 

0.83843 

0.94355 

1.15717 

0.0 

0.01443 

0.03093 

0.04075 

3.50 

0.90399 

0.91639 

0.97178 

1-.09111 

0.0 

0.01275 

0.02786 

0.03839 

4.00 

0.94397 

0.95081 

0.98222 

1.05509 

0.0 

0.01163 

0.02601 

0.03777 

5.00 

0.97560 

0.97815 

0.99037 

1.02238 

0.0 

0.01036 

0.02400 

0.03799 

6.00 

0.98687 

0.98803 

0.99373 

1.00994 

0.0 

0.00971 

0.02298 

0,. 03860 

8.00 

0.99455 

0.99490 

0.99664 

1.00208 

0.0 

0.00909 

0.02201 

0.03962 

10.00 

0.99705 

0.99719 

0.99789 

1.00018 

0.0 

0.00881 

0.02158 

0.04025 

15.00 

0.99891 

0.99894 

0.99908 

0.99954 

0.0 

0.00854 

0.02116 

0.04099 







Anticipating the results of numerical experiments, to be discussed, the 
phase-accurate optimum value of v”^ = introduces entirely too much arti- 
ficial diffusion, for both linear and non-linear example equation systems. 

In addition, these linearized analyses are valid only for the algorithm using 
linear interpolation, and performance assessment with at least quadratics is 
required. As an indication of expected performance, the quadratic interpola- 
tion basis yields formally fourth order accurate difference representation 
for both convection and diffusion differential operators at the elemental 
vertex nodes of a uniform discretization of x“. The additional required as- 
sessment of the tensor product algorithm basis is also facilitated by numeri- 
cal experiment. 


NUMERICAL RESULTS 

One-Dimensional Solutions 

The primary requirement is to assess acceptable bounds on v > 0 that 
facilitate accurate solutions without introduction of excessive artificial 
diffusion. An appropriate example for examining the important non-linearity 
of the momentum equation (2) is the inviscid form of Burgers equation 

L(u)=||-+u||=0 (42) 

For an initial condition corresponding to a square wave, see Figure 2a), the 
exact solution to equation (42) is propagation of the original wave form 
parallel to u with a celerity of igu. Figures 2b)-c) show the results obtained 
from the dissipative finite element algorithm for k = 1, C = 0.125 and 
= 1 ^. The propagation speed of the wave is exactly correct, ie., the 





wave numerical distribution is identically repeated every 16 time steps, ie. 
2/C, The solution approximation to the original square wave is excellent, 
ie., the step remains interpolated across one element domain only, the 
leading phase error is only 0.35J, while the lagging phase error extremum is 
1.3%. For this value of v, there is no perceptible diffusion of the step; 
however,' diffusion is introduced when the "optimal" linear analysis value of 
V ^ = vTF is used, and the square wave becomes interpolated across three 
element domains. 

Figure 3 summarizes the influence of Courant number (time intergration 
step size) and level of v on the square wave solution. The results in Figure 
3a) were obtained for C = 0,5 and v"^ = /SU". The fidelity of ' the original 
square wave is excellently maintained with no evidence of numerical' diffusion. 
The leading phase error is reduced to 0.1%, compare to Figure 2b), while the 
lagging error extremum is increased to 1.9%. The algorithm is stable to unit 
Courant Number; the fidelity of the original square wave is degraded further 
(lagging error extremum is 4.6%), as induced by the truncation error associ- 
ated with the trapezoidal rule. However, the wave remains interpolated 
across only one finite element. For comparison, setting v~^ = and C = 1.0 
yields the results shown in Figure 3c). The lagging error peaks have become 
completely diffused, the leading error peak is 1.9%, and the wave has be- 
come interpolated over four element domains. 

These results verify that the developed dissipative implicit finite 
element algorithm, employing linear interpolation bases, exhibits excellent 
accuracy control for the sample non-linear problem. Additional numerical 
tests have quantized the relative importance of the derived matrix structures. 
By and large, all modification degrade performance of the basic algorithm. 
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c) v‘^ = C = 1.0 

F1g. 3. Solution of l-D Burger's Equation, Various Courant 

Numbers and Viscosity, Linear Elements, nAt = 100 
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For example, the results in Figure 4a) were obtained following replacement 
of the diffusion term v{U*^}g v;ith the scalar v{l}^ in equation (30). The 
results are clearly inferior to the comparison test. Figure. 2. For Figure 
4b), this modification was retained, and in addition, the finite element- 
derived initial value matrix structure, equation (49), was collapsed, to the 
finite difference diagonal form. A large lagging dispersion error peak is 
introduced, and these results are definitely poorer than the comparison case. 
Figure 3a). Figure 4c) corresponds to the exact duplicate of the Beam and 
Warming (ref. 13)' implicit finite difference algorithm with an added fourth- 
order dissipation term, ie., a fourth order accurate finite difference equi- 
valent of a viscosity term. The value of y, the artificial diffusion coef- 
ficient, was numerically optimized, and Figure 4c compares almost exactly 
with Figure Ic of reference 13. By comparison, these results are clearly much 
poorer than those of Figure 3. Hence, these tests firmly quantize the super- 
ior performance of the finite element-based algorithm for this test case. 

Similar results are obtained for the quadratic element embodiment of 
the algorithm, as obtained setting k = 2 in equation (20), Figure 5a) illus- 
trates the square wave after 80 time steps at C = 0.125 and for v”^ = ✓T5'. 
Since ‘the quadratic possesses a non-vertex node, the original wave (inter- 
polated across one element) possesses a nodal mid-value. The quadratic al- 
gorithm maintains an adequately accurate representation of the original wave, 
see Figure 5a), with lagging and leading dispersion error extrema of 4%. 

The phase celerity is again exactly correct, as the numerical solution is 
repeared every 16 time steps. For Figure 5b), the dissipation term was 
seal ari zed, as discussed for the k = 1 solution in Figure 4b), which induces 
sufficient additional dissipation to eliminate the dispersion error and 
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a) k = 1; S = 0; {1}; C = .125 



b) k = 1; S = 2; v"^= v/so {!}; C = 0.5 
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Solution of 1-D Burger's Equation, Various Assembly and 
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Fig. 5, Solution of 1-D Burger's Equation, Quadratic Elements, • 
Various Assembly and Viscosity, nAt - 80 
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smear the wave over two element domains. The additional diagonalizing of 
the initial value matrix does not markedly alter this solution, Figure 5c, 
in contrast to the linear element solution results. 

The influence of the dissipative mechanism within the derived algor- 
ithm is less demonstrative for a one-dimensional linear equation solution, 
eg. the continuity equation (1). A test case is advection of a cosine wave 
by a constant imposed velocity, for which the theoretical analysis is exact. 
Figure 6a) shows the initial-condition, and Fig. 6b) the non-dissipative 
linear (v = 0, k = 1) finite element solution following propagation over 
three wave-lengths. The dashed curve is the exact solution, and for C = 0.4 
the finite element result agrees almost exactly. A modest leading, phase 
error is exhibited, a trailing 2-3Ax wave is induced by the relatively 
poorer phase accuracy, and the peak value remains at 100%. The diagonal 
initial-value matrix form yields. exactly the Crank-Nicolson algorithm. Fig- 
ure 6c) shows the corresponding results, which are substantially poorer in 
comparison. As predicted by the theoretical analysis, phase fidelity is much 
poorer, to the extent that the marginally diffused peak celerity becomes in 
substantial error. In contrast to the non-linear square wave test case, 
this solution if continued will produce totally erroneous trash. 

Table 2 summarizes the accuracy of the non-dissipative linear element 
algorithm solutions as a function of integration time-step, ie. Courant num- 
ber. A peak value of 100 is retained up to C 0.5, and solution dissymetry 
progressively increases with larger Courant number. The last column denotes 
the magnitude of the largest trailing wave peak, which always occurs imme- 
diately behind the cosine wave. (Not until C = 1.4 is the finite element 
solution distribution similar in appearance to the Crank-Nicolson result for 
C = 0.4.) Table 3 summarizes the influence of the dissipation factor level 
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Table 2 


Influence of Courant Number on 
Cosine Hill Distribution, k = 1, v = 0 


Courant 

No, 

Solution Distribution 

Dispersion 

Error 

Peak 

Analytical 

15 

50 

85 

100 

85 

50 

15 

0 

0.1 

16 

50 

85 

100 

84 

47 

15 

-2 

■ 0.2 

16. 

50 

85 

100 . 

83 

46 

15 

-2 

0.4 

17 

52 

87 

100 

81 

45 

15 

-3 

0.5 

18 

54 

90 

lOp 

78 

43 

16 

-5 

0.6 

20 

56 

90 

98 

75 

42 

16 

-6 

0.7 

21 

59 

91 

96 

73 

41 

17 

-7 

0.8 

23 

61 

92 

94 

70 • 

40 

17 

-8 
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Table 3 


Influence of Dissipation Level on Solution 
For Cosine Hill, C = 0.4, k = 1 


Dissipation 
Level (v) 



Sol ution 

Distribution 



Dispersion 

Error 

Peak 

Analyti cal 

15 

50 

85 

100 

85 

50 

15 

0 

0.0 

17 

52 

87 

100 

81 

45 

15 

-3 

0.1 E-03 

18 

52 

85- 

98 

80 

45 

17 

-3 

0.1 E-02 

19 

52 

85 

98 

80 

45 

17 . 

-3 

0.1 E-01 

19 

52 

84 

97 

79 

45 

17 

-3 

0.2 E-01 

20 

52 

84 

95 

78 

46 

18 

-2 

0.5 E-01 

21 

52 

81 

92 

76 

46 

20 

-1 

0.9 E-01 

23 

52 

78 

87 

74 

47 

22 

0 

0.1291 

25 

51 

75 

84 

72 

47 

23 

0 
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V on phase accuracy for the test case at C = 0.4. Levels of v < 0.01 do 
not measureably alter the solution. Increasing v < 0.01 tends to pro- 
gressively symmetrize the solution while continuously adding diffusion of 
the peak level. At v = 0.1291, which is one-half the linear analysis 
optimum value, the solution is nearly symmetric, but the level of peak 
diffusion is unacceptably large. Therefore, it is confirmed that the 
sixth-order phase accuracy is unacceptable in terms of induced artificial 
diffusion. 

The corresponding results- for the quadratic element form of the 
algorithm are summarized in Tables 4 and 5. Up to C = 0.5, the non- 
dissipative algorithm produces essentially identical results using k = 1 
and 2, compare Tables 2 and 4. For larger Courant number, the inaccuracy 
and phase shift of the quadratic form becomes progressively poorer in 
comparison. Comparing Tables 3 and 5, there is little performance 
difference between k = 1 and k = 2 for v > 0 at C = 0.4 for this test case. 
Definitive. differences will result for multidimensional solutions, however. 
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Table 4 


Influence of Courant Number on 



Cosine 

Hill 

Distribution, 

k = 2, 

V = 0 



Courant 

No. 


Solution 

Distri bution 



Dispersion 

Error 

Peak 

Analytical 

15 

50 

85 

100 

85 

50 

15 

0 

0,1 

14- 

49 

85 

98 

87 

50 

14 

-1 

• 0.2 

14 

50 

86 

99 

86 

. 49 

14 

-1 

0.3 

15 

52 

87 

99 

84 

47 

13 

-2 

0.4 

16 

53 

89 

99 

82 

45 

13 

-3 

0.5 

19 

58 

92 

98 

77 

41 

14 

-7 

0.6 

21 

60 

94 

97 

74 

40 

15 

-9 

0.7 

26 

66 

96 

91 

67 

38 

16 

-11 

0.8 

26 

67 

98 

92 

67 

37 

16 

-13 
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Table 5 


Influence of Dissipation Level on Solution 
For Cosine Hill, C - 0.4, k = 2 


Dissipation 
Level (v) 


Solution 

Distribution 



Dispersion 

Error 

Peak 


15 

50 

85 

100 

85 

50 

15 

0 

o.p 

16 

53 

89 

99 

82 

45 

13 

-3 

0:i E-3 

16 

53 

89 

98 

82 

45 

13 

-3 

0.1 E-2 

17 

53 

88’ 

98 

82 

45 

14 

-3 

0.1 E-1 


•53 

85 

94 

79 

46 

16 

-2 

0.2 E-1 

21 

53 

82 

90 

77 

46 

19 

-1 

0.5 E-1 

26 

51 

73 

79 

70 

47 

24 

0 

0.9 E-1 

30 

49 

65 

70 

63 

46 

28 

0 ' 

0.1291 

31 

47 

59 

63 

58 

45 

30 

0 
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Multi -Dimensional Solutions 


The one-dimensional solutions have quantized acceptable levels of 
V > 0 for which adequate accuracy can be maintained within acceptable dif- 
fusion levels. The present requirement is embedding the tensor-product 
basis into a multi -dimensional algorithm and evaluate non-linear and linear 
solutions. Figure 7 illustrates, for two-dimensions, the manner in which 
fractional steps are employed for the grid sweeps associated with use of 
the tensor product basis. For linear interpolation. Figure 7a), the ele- 
ments of the linearly independent. cardinal basis are contracted on the first 
sweep (parallel to axis) with node numbers 


{1, 2; 2. 3; 3, 4; 29, 30; 31, 32; 32, 33; ...} 


to form the vectors on the interpolation subspace. Within the fractional 
steps concept, the first (inner) tensor matrix statement in equation (25) is 
iterated to convergence of the second sweep, parallel to the 

X 2 axis, the elemental contraction vectors' are ordered 
{1, 31; 31, 61; ...; 200, 201; 2, 32; ...} 
and the independent variables in the interpolation basis are linearly depen- 
dent upon X 2 . The iteration of the interior matrix product then yields the 
converged value of The extension to the third direction is obvious, 

whereupon the matrix iteration converges to the solution at time-step 


The linear element tensor product algorithm is rather comparable to a 
finite difference alternating direction framework. The quadratic basis al- 
gorithm is somewhat more complex. Referring to Figure 7b), for the first 
sweep the contraction nodal vector is 

{1, 2, 3; 3, 4, 5; 28, 29, 30; 31, 32, 33; ...} 
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b) Quadradic Basis 


Fig. 7. Illustrations for Tensor Product Basis Finite Element Algorithm 







The even numbered nodes 2 , 4, 6, 32 , ... are non-vertex in the 

first sweep, which yields the solution the second sweep 

parallel to X 2 , the contraction vector becomes 

{1, 31, 61; 61, 91, 121; ...; 2, 32, 62; ...} 

The odd-number nodes 31, 91, 33, etc, are now non-vertex, while 2, 62, 

4, etc. have become vertex, and the solution is Mote that only 

nodes 32, 34, 92, etc. are always non-vertex for a two-dimensional prob- 
lem. Hence, the tensor product algorithm employs an interpolation basis 
essentially comparable to the non-Serendipity multi -dimensional quadratic 
with an interior (centroidal) node. The extension of the quadratic ten- 
sor product algorithm to three-dimensions is again direct. To avoid the 
one-sided accumulation of round-off error, the tensor matrix algorithm 
sweeps are sequenced with each coordinate direction cyclically alternating 
position within equation (25). 

The comparable non-linear test case is the two-dimensional inviscid 
Burgers equation system 

= ■(«) 

Figure 8a shows the initial condition for a square wave impinging on the 
upper left corner of the solution domain. The initial step distribution is 
interpolated across one element domain only; the correct solution is pure 
advection of the initial distribution, parallel to thedomain diagonal with 
celerity and maintenance of unit and zero level plateaus. Vanishing 

normal derivative boundary conditions were applied everywhere except for 
those nodes possessing non-zero initial values which were fixed. Figures 
8b) -c) show the linear element computed results for C = 0.125 and v"^ = v'SF. 
They are excellent approximations to the correct solution away from the 
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c) Downstream, nAt = 120 


Fig. 8. Solution of Two-Dimensional Burgers Equation, Linear Tensor 
Product Finite Element Algorithm, C = .125, 
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boundaries. Thereupon, the numerical solution generates an acceptable 
approximation to the vanishing normal derivative (which cannot be explicitly 
enforced, see equation (6), since Jl(q) = 0 is a homogeniety in equation 
{9)). The center wave celerity is exact, and the dispersion error-induced 
peaks are nominally + 4°^ and are limited to the immediate vicinity of the 
front. For comparison. Figure 9 shows the identical test case but with 

V s 0 within the algorithm. The importance of the dissipative mechanism 
is graphically evident. 

As was done for the one-dimensional cases, comparison results are ob- 
tained with finite difference matrix structure modifications. The di ag- 
onal izati on of the Initial -value matrix degrades the algorithm to second- 
order accuracy. The concomitant artificial diffusion is further enhanced 
by reducing the phase selectivity of the dissipation action by setting 

V {U} ^ V {1}. The resultant solution. Figure 10 has diffused the front 
over approximately six element domains. The wave front is considerably 
sharpened, see Figure 11, by restoring the theoretical matrix statement 
v{U}, with an associated increase (to \ 1 %) of the lagging phase error peak. 
These multi -dimensional results appear quite comparable to the finite dif- 
ference experience discussed in Figure 4. In all cases, the restructuring . 
of the basic theoretical statement degrades solution performance. 

Figure 12 shows the corresponding test case solutions generated by 
the quadratic tensor product algorithm for C = 0.125 and v"^ = /SO". The 
solutions are basically identical to those obtained by the linear element 
form, see Figure 8, with only modest leading and lagging error peaks, ex- 
cellent plateaus and good approximations to vanishing normal boundary 
gradients. In particular, the one-dimensional wave form in Figure 5a) is 
almost identical with the intersection of the two-dimensional solution with 
the plane y = 0, Figure 12b), 
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a) Downstream, nAt = 60 



Fig. 9. Solution of Two-Dimensional Burgers Equation, Linear Tensor 
Product Finite Element, C = .125, v s 0 
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a) Downstream, nAt = 60 



b) Downstream, nAt = 120 


Fig. 10. Solution of Two-Dimensional Burgers Equation, Diagonalized Linear 
Tensor Product Algorithm, C = 0.125, v"^ = /SO" ({1}), 
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a) Downstream, nAt = 60 



b) Downstream, nAt =120 


Fig. 11. Solution of Two-Dimensional Burgers Equation, Diagonalized Linear 
Tensor Product Algorithm, C = 0.125, v ^ {{U}). 
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Solution improvement with the quadratic tensor product basis occurs 
for the multi -dimensional solution of the linear equation. The two-dimen- 
sional test case equivalent is the "cosine hill", obtained by rotating the 

J 

one-dimensional cosine v/ave about the peak. A demanding test is pure ad- 
vection by an imposed solid body rotation, cf. reference 16. The correct 
solution is exact propagation of the initial distribution. Fibures 13a) -c) 
illustrate the linear element solution results for C = 0.15 and v e 0. After 
one-full turn, the peak value would remain at 100 if the phase distortion did 
not move it off a node location, and the lagging phase error peak is 10%. 

Figure 14 quantizes the data; Figure 14a) shows the initial condition distri- 
bution, which is identical with the exact solution following one full revolu- 
tion. Figure 14b) shows the k = 1 results obtained after one rotation at 
C = 0.15, and the circled value corresponds to the correct peak location. 

The generally lagging phase has retarded the computed peak about one-half cell, 
and the phase dispersion error is firmly quantized. Figure 15a) shows the 
comparison solution, obtained with the non-dissipative k = 2 algorithm form 
at C = 0.15, the accuracy of which is excellent. The solution distortion due 
to lagging phase is nominally absent, the lagging wake, peak is -2%, and the 
solution peak is undiffused (it is actually modestly increased) and occurs 
at the exact nodal location. Doubling the integration step size to C = 0.3 
produces the solution shown in Fig. 15b), which is essentially comparable to . 
the linear element solution at C = 0.15. Due to the pentadi agonal Jacobian, 
the quadratic algorithm is about 16% slower than the linear; therefore, the 
net CPU savings for the k = 2 solution at double the Courant Number is about 
35%. The linear element solution at C = 0,3 further decreased the peak value 
to 82, from the 93 in Figure 14b), and the corresponding enhancement of dis- 
persion error produced an unacceptable -17 in the trailing wake. 
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Fig, 13. Advection of Cosine Hill In Solid-Body Rotation Velocity 
Field, C = 0.15, k = 1. 
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a) Initial Condition, Exact Final Solution 
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b) Linear Tensor Product Solution, One Rotation, C = 0.15, v = 0. ‘ 
Fig. 14. Advection of Cosine Hill in Solid Body Rotation 
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Fig. 15. Advection of Cosine Hill in Solid Body Rotation 
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A three-dimensional solution for which an exact solution is known is 
pure advection of the "cosine-sphere", the three-dimensional equivalent of 
the "cosine-hill." The selected case is linear convection along the domain 
diagonal, and Figure 16 visualizes the solution in terms of particle density 
distributions. Here, the lagging phase of the k = 1 algorithm progressively 
sheds particles into the distribution wake. The accuracy of the tensor pro- 
duct algorithm is quantized in the remaining figures, which are printouts of 
the computed two dimensional distributions in the three Xs-planes centered 
about the exact solution mid-plane. For comparison. Figure 17 contains the 
initial distributions; the correct solution is pure translation to the upper 
right corner, and preservation of all symmetries. 

The result for the non-dissipative linear algorithm, as obtained for 
C = 2/3, are shown in’ Figure 18. The peak level is almost retained (99), 
and the lagging phase distortion produces overall lower solution levels in 
the upper plane. Fig. 18c) compared to the lower plane. The dispersion error 
produces the evidenced trailing wakes with peaks of - 1 %. The solution alter- 
ations produced by introducing v > 0 are summarized in Figure 19, which are 
printouts of central-plane distributions for v = 0.006 and 0.012. The dis- 
persion error peaks are modestly reduced, with the corresponding decrease in 
peak level. 

The comparison results for the quadratic algorithm at C = 2/3 are shown 
in Figure 20 for the non-dissipative form. . The peak value is enhanced, as 
occurred for the two-dimensional solution, and overall symmetries are con- 
siderably improved over the linear element results. The dispersion error peaks 
in the trailing wake are also modestly higher in comparison. In Figure 21a), 
setting v = 0.06 diffuses the peak to 96%, and reduces the wake error by 
nominally half while retaining the symmetry preservation of the v = 0 solution. 
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Fig. 17. Three-Dimensional Planar Initial Distributions at Center 
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Fig. 20. Computed Centroidal Plane Distributions for Three-Dimensional- 
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Fig. 21. Computed Centroidal Plane Distributions for Three-Dimensional 
Solutions, k = 2, C = 2/3. 


- 50 - 



Increasing v = 0.012 diffuses the peak further to 87% while further improving 
the wake error, Figure 21. b). AH the three-dimensional results further con- 
firm that acceptable levels of dissipation v are appreciably below the linear 
analysis optimum value. 

The computer requirements for the three-dimensional solutions encompass 
those of all the lower dimensional cases. A total of 120,000 words of core 
were sufficient to execute both the linear and quadratic algorfthm. forms on 
17^ mesh. The Jacobian LU decomposition, for the quadratic requires about 
tvricG the storage of the linear, but the Jacobian core constitutes less than 
1 % of total storage using the tensor product algorithm. The CPU to execute 
one sv/eep of the quadratic algorithm form is approximately 15 - 20% larger than 
to execute the linear element sweep. The linear element three-dimensional test 
case requires less than 1 minute of CPU on an IBM 360/195 computer for execu- 
tion. 

SUMMARY AND. CONCLUSIONS 

An accurate and efficient tensor product basis finite element solution 
algorithm is established for application to convection dominated flow field 
predictions. The intrinsic fourth-order accuracy is enhanced using a dissi- 
pative formulation to modify phase error-induced oscillations and instabili- 
ties, Embedding the formulation within the method of fractional steps yields 
a core-efficient procedure for implicit integration of the resultant ordinary 
differential equation systems. 

The results of numerical experiments, for linear and non-linear model 
partial differential equations, firmly quantize the performance and accuracy 
of the developed algorithm. In particular, practically acceptable levels of 
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numerical diffusion demand the dissipation parameter be selected much smaller 
than an "optimum" linear analysis evaluation. The overall performance of 
the quadratic tensor product basis algorithm is modestly superior to the 
linear basis form, although the latter functions quite well for an exceed- 
ingly simple formulation. Minor modifications to the derived matrix structures 
for the linear algorithm yields familiar finite difference forms, the perform- 
ance of which appears inferior based upon the results presented. The deve- 
loped finite element algorithm should find wide application in computational 
fluid dynamics. 
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